PDE1B, a potential biomarker associated with tumor microenvironment and clinical prognostic significance in osteosarcoma

PDE1B had been found to be involved in various diseases, including tumors and non-tumors. However, little was known about the definite role of PDE1B in osteosarcoma. Therefore, we mined public data on osteosarcoma to reveal the prognostic values and immunological roles of the PDE1B gene. Three osteosarcoma-related datasets from online websites were utilized for further data analysis. R 4.3.2 software was utilized to conduct difference analysis, prognostic analysis, gene set enrichment analysis (GSEA), nomogram construction, and immunological evaluations, respectively. Experimental verification of the PDE1B gene in osteosarcoma was conducted by qRT-PCR and western blot, based on the manufacturer's instructions. The PDE1B gene was discovered to be lowly expressed in osteosarcoma, and its low expression was associated with poor OS (all P < 0.05). Experimental verifications by qRT-PCR and western blot results remained consistent (all P < 0.05). Univariate and multivariate Cox regression analyses indicated that the PDE1B gene had independent abilities in predicting OS in the TARGET osteosarcoma dataset (both P < 0.05). GSEA revealed that PDE1B was markedly linked to the calcium, cell cycle, chemokine, JAK STAT, and VEGF pathways. Moreover, PDE1B was found to be markedly associated with immunity (all P < 0.05), and the TIDE algorithm further shed light on that patients with high-PDE1B expression would have a better immune response to immunotherapies than those with low-PDE1B expression, suggesting that the PDE1B gene could prevent immune escape from osteosarcoma. The PDE1B gene was found to be a tumor suppressor gene in osteosarcoma, and its high expression was related to a better OS prognosis, suppressing immune escape from osteosarcoma.

www.nature.com/scientificreports/Phosphodiesterase 1B (PDE1B) belonged to the phosphodiesterase (PDE) family and PDE1 subfamily, including three main PDE1 isoforms of PDE1A-1C and being stimulated by a calcium-calmodulin complex 10,11 .Zang et al. shed light on that PDE1B had the ability to regulate exosome biogenesis and autophagic flux in microglia for neuroprotection under ischemic stroke conditions 12 .McQuown et al. revealed that PDE1B could serve as a negative regulator of memory in the hippocampus, making it a promising target for memory impairment 13 .In clear cell renal cell carcinoma, PDE1B was found to be the target of miR-5701 in promoting the tumor cells' apoptosis 14 .Chen et al. established a novel tumor microenvironment (TME)-related signature, including PDE1B, and this signature could predict survival prognosis and therapeutic responses for colon cancer patients 15 .In osteosarcoma, only two gene signatures, including PDE1B, were available, and these two signatures were related to osteosarcoma patients' overall survival (OS) or metastasis 16,17 .Currently, little was known about the definite role of PDE1B in osteosarcoma.Therefore, the present article was aimed at mining public data on osteosarcoma to reveal the prognostic values and immunological roles of the PDE1B gene in osteosarcoma, with the assistance of experimental verifications.Our results were anticipated to provide operational targets for future clinically personalized treatment targets in osteosarcoma.

Gene screening and data processing
A total of three osteosarcoma-related datasets were utilized for further analysis, containing the GSE28424 osteosarcoma cell line dataset (https:// www.ncbi.nlm.nih.gov/ geo/ query/ acc.cgi? acc= GSE28 424), including 19 osteosarcoma cells and 4 control cells; the GSE33382 osteosarcoma tissue dataset (https:// www.ncbi.nlm.nih.gov/ geo/ query/ acc.cgi? acc= GSE33 382), including 84 osteosarcoma tissues and 3 control tissues; and the Therapeutically Applicable Research to Generate Effective Treatments (TARGET) osteosarcoma tissue dataset (https:// www.cancer.gov/ ccg/ resea rch/ genome-seque ncing/ target), including 86 osteosarcoma tissues.With the help of the R "limma" version 3.58.1 package, differential expression analyses were obtained from the GSE28424 and GSE33382 datasets with cut-off values of P values below 0.05 for identifying differently expressed genes via the Wilcoxon signed-rank test.By using the R "survival" version 3.5-7 package, survival differences were obtained from the TARGET osteosarcoma dataset with cut-off values of P values below 0.05 for identifying OS-significant genes.A Venn diagram was conducted to obtain intersecting genes from these differently expressed and OSsignificant genes.Finally, a total of 10 intersected genes were obtained, including the PDE1B gene.Based on our experimental results, the PDE1B gene was selected, and the TARGET osteosarcoma dataset was further utilized for analyzing the roles of the PDE1B gene in osteosarcoma.K-M survival analysis was conducted based on the medium expression of the PDE1B gene, which classified high-and low-risk groups.These analyses were all conducted by the R 4.3.2software.

Experimental verifications of PDE1B gene in osteosarcoma
Based on the manufacturer's instructions, quantitative real-time PCR (qRT-PCR) and western blot were conducted to verify the expression levels of the PDE1B gene in osteosarcoma cell lines (MG-63, SW1353, 143B, and hFOB) purchased from Procell (Wuhan, China) using the methods previously described 18,19 .Therein, qRT-PCR was utilized to verify the mRNA expression levels of the PDE1B gene in osteosarcoma cell lines, and the primers were displayed as follows: human-PDE1B-F: 5ʹ-GAG GCT CCA TCC GAC CAA T-3ʹ, human-PDE1B-R: 5ʹ-GCA CCC TTG ACC CTA ACC C-3ʹ; human-ACTB-F: 5ʹ-TAG TTG CGT TAC ACC CTT TCTTG-3ʹ, human-ACTB-R 5'-TGC TGT CAC CTT CAC CGT TC-3ʹ.Western blot was utilized to verify the protein expression levels of the PDE1B gene in osteosarcoma cell lines, and the antibodies were obtained from Abcam (ab182565; https:// www.abcam.cn/).

Associations between PDE1B gene and clinical factors in osteosarcoma
As previously described 20,21 , difference analysis, Cox regression analysis, and nomogram were utilized to reveal the associations between the PDE1B gene and clinical factors in osteosarcoma, respectively.Therein, boxplots were applied to display the PDE1B gene expression distribution in different clinical factors (gender of female and male; race of African, Asian, and White; first event of relapse, censored or none; disease at diagnosis of metastatic or nonmetastatic; primary tumor site of arm/hand/pelvis, and leg/foot).Moreover, it is statistically called censored if the outcome event does not occur at the specified end time due to loss of visit, death, failure to recover, etc. Univariate and multivariate Cox regression analyses were utilized to find clinical factors with independent abilities in predicting OS for osteosarcoma.A nomogram was also put in place to intuitively forecast the survival probabilities of osteosarcoma patients.C-index, ROC curves, and calibration plots were utilized to assess the performance of the established nomogram in osteosarcoma.These analyses were all conducted by the R 4.3.2software.

Pathway predictions and immunological associations
To obtain PDE1B related pathways in osteosarcoma, gene set enrichment analysis (GSEA) was applied in high-PDE1B and low-PDE1B subgroups by the GSEA 4.0.0 software with default parameters as previously described 22,23 .To reveal the associations between PDE1B gene expression and immunity in osteosarcoma, tumor immune cells infiltration levels and tumor microenvironment were assessed based on the medium expression of the PDE1B gene, classified into high-and low-risk groups.Tumor immune cells infiltration levels were calculated by the "CIBERSORT" version 1.03 R script to assess the infiltration levels of 22 kinds of immune cells in high-PDE1B and low-PDE1B subgroups 24,25 .Tumor microenvironment was calculated by the "ESTIMATE" algorithm of the R "estimate" version 1.0.13package to assess immune, ESTIMATE, and stromal scores in high-PDE1B and low-PDE1B subgroups 26,27 .The Tumor Immune Dysfunction and Exclusion website (TIDE; http:// tide.dfci.harva

Ethics approval and consent to participate
This study was approved by the Institutional Research Ethics Committees of Affiliated Hospital of Nantong University.

The expression of PDE1B in osteosarcoma
The whole analyzing process was detailed in Fig. S1.The Venn diagram showed that a total of 10 intersected genes were obtained from three osteosarcoma-related datasets (GSE28424, GSE33382, and TARGET osteosarcoma datasets), including the PDE1B gene (Fig. 1A).Therein, PDE1B was differently expressed in the GSE28424 and GSE33382 datasets, with the cut-off values of the P value below 0.05 (Fig. 1B,C).Moreover, PDE1B was significantly associated with OS in the TARGET osteosarcoma dataset, with a cut-off value of P value below 0.05 (Fig. 1D).Experimental verifications by qRT-PCR and western blot results indicated that both PDE1B gene mRNA and protein expression levels were lowly expressed in osteosarcoma cell lines (all P < 0.05; Fig. 1E-G).The above-mentioned results showed that PDE1B was lowly expressed in osteosarcoma and significantly associated with OS.

Clinical factors' relationships with PDE1B gene in osteosarcoma
Figure 2 presented the PDE1B gene expression distribution in different clinical factors (gender, race, first event, disease at diagnosis, and primary tumor site).Our results indicated that the PDE1B gene was remarkably related to metastasis (P < 0.001; Fig. 2D) in the TARGET osteosarcoma dataset, while it was not associated with gender, race, first event, or primary tumor site (all P > 0.05).Univariate and multivariate Cox regression analyses indicated that the PDE1B gene had independent abilities in predicting OS in the TARGET osteosarcoma dataset (both P < 0.05; Fig. 3 and Table 1).To further forecast the 1-, 3-, 5-year OS probabilities of osteosarcoma patients intuitively, a nomogram was established based on several clinical factors (gender, race, first event, disease at diagnosis, and PDE1B gene expression) (Fig. 4A).C-index and 1-, 3-, 5-year AUCs of the PDE1B-based nomogram were 0.867, 0.892, 0.899, and 0.873 in the Target osteosarcoma dataset, respectively (Table 2).1-, 3-, 5-year calibration plots indicated the consistency between actual and ideal results, suggesting a good performance of our established nomogram in osteosarcoma (Fig. 4B-D).The above-mentioned results showed that PDE1B was significantly associated with metastasis and had independent abilities in predicting OS in osteosarcoma.

PDE1B related pathways and its associations with immune checkpoints genes, m6A genes, VEGF pathway genes in osteosarcoma
GSEA was applied by us to reveal PDE1B related pathways in osteosarcoma, and our outcomes showed that PDE1B was markedly linked to calcium, cell cycle, chemokine, JAK STAT, and VEGF pathways (Fig. 5 and Table 3).We further explored the immune checkpoint gene, m6A gene, and VEGF pathway gene expression levels in low-and high-PDE1B subgroups (Fig. 6).Our results reported that the PDE1B gene was significantly   Table 1.Prediction of overall survival by using univariate and multivariate analyses of PDE1B expression level and clinicopathological variables in Target osteosarcoma dataset.Bold font means P < 0.05.HR hazard ratio, HR.95L lower limit of HR 95% confidence interval, HR.95H upper limit of HR 95% confidence interval, Gender female or male, Race African, Asian, or White; Dis_diagnosis metastatic or nonmetastatic, Tumor_site arm/hand/pelvis or leg/foot.involved with immune checkpoint genes and VEGF pathway genes.For m6A regulator genes, only RBMX and RBM15 were firmly related to PDE1B gene expression.The above-mentioned results showed that PDE1B was significantly associated with five pathways, immune checkpoint genes and VEGF pathway genes in osteosarcoma.

Associations between PDE1B gene expression and immunity in osteosarcoma
To reveal the associations between PDE1B gene expression and immunity in osteosarcoma, tumor immune cells infiltration levels and tumor microenvironment were assessed based on the medium expression of the PDE1B gene, classified into high-and low-risk groups.Tumor immune cells infiltration levels were calculated by the CIBERSORT algorithm to assess the infiltration levels of 22 kinds of immune cells in high-PDE1B and low-PDE1B subgroups (Fig. 7A).The infiltration levels of T cells follicular helper, T cells gamma delta, Macrophages M0, and Neutrophils were differently expressed in high-PDE1B and low-PDE1B subgroups (all P < 0.05; Fig. 7B).
Further correlation analyses showed that PDE1B gene expression was significantly associated with T cells gamma delta or Macrophages M0 cells infiltration levels (all P < 0.05; Fig. 7C,D).The tumor microenvironment was calculated by the ESTIMATE algorithm to assess immune, ESTIMATE, and stromal scores in high-PDE1B and low-PDE1B subgroups.Our results presented that immune, ESTIMATE, and stromal scores were differently expressed in high-PDE1B and low-PDE1B subgroups (all P < 0.05; Fig. 7E).Moreover, correlation analyses showed that PDE1B gene expression was significantly associated with immune, ESTIMATE, and stromal scores (all P < 0.05; Fig. 7F-H).The above-mentioned results showed that PDE1B was significantly associated with immunity in osteosarcoma.

Evaluating osteosarcoma patients' PDE1B gene expression for immunotherapies by TIDE dataset.
The TIDE algorithm was also applied to forecast immune responses by evaluating patients' PDE1B gene expression for immunotherapies.The outcomes of us shed light on that patients with high-PDE1B expression would  www.nature.com/scientificreports/have a lower TIDE score and a lower T cell dysfunction score (Fig. 8).The above-mentioned results showed that patients with high-PDE1B expression would have a better immune response to immunotherapies than those with low-PDE1B expression in osteosarcoma.

Discussion
Osteosarcoma was characterized by poor prognoses, frequent recurrences, as well as distant metastasis properties, seriously affecting children and adolescents aged 15-19 years old 30 .Despite the progress in treatments for osteosarcoma containing surgery, adjuvant chemotherapy, radiotherapy, drug therapy, immunotherapy, and so on, more than one in five osteosarcoma patients were diagnosed with metastases, leading to a five-year survival rate of less than 30% 31 .The PDE1B gene belonged to the PDE family, and this gene had been found to be involved in various diseases, including tumors and non-tumors.Currently, little was known about the definite role of PDE1B in osteosarcoma.Therefore, we mined public data on osteosarcoma to reveal the prognostic values and immunological roles of the PDE1B gene, providing operational targets for further clinically personalized treatment targets in osteosarcoma.
In this article, we firstly mined three osteosarcoma related datasets (GSE28424, GSE33382 and TARGET osteosarcoma datasets) and identified a total of 10 intersected genes by Venn diagram, including the PDE1B gene.As reported, Tan et al. identified and confirmed the roles of ABCG8, PDE1B, and LOXL4 in osteosarcoma for predicting OS 16 .Wen et al. also found a three-gene signature, including MYOM2, PDE1B, and COCH genes,  www.nature.com/scientificreports/ that could predict the OS of osteosarcoma 17 .In our article, the PDE1B gene was discovered to be lowly expressed in osteosarcoma, and its low expression was significantly associated with poor OS.Experimental verifications by qRT-PCR and western blot results remained consistent.In the associations with clinical factors, the PDE1B gene was remarkably related to metastasis.Univariate and multivariate Cox regression analyses indicated that the PDE1B gene had independent abilities in predicting OS in the TARGET osteosarcoma dataset.We also established a nomogram based on several clinical factors (gender, race, first event, disease at diagnosis, and PDE1B gene expression) to further forecast the 1-, 3-, 5-year OS probabilities of osteosarcoma patients intuitively with good performance.GSEA revealed that PDE1B was markedly linked to the calcium, cell cycle, chemokine, JAK STAT, and VEGF pathways.Moreover, PDE1B was found to be markedly associated with immunity, including the T cells gamma delta, Macrophages M0 cell infiltration levels, and the tumor microenvironment.Further TIDE algorithms shed light on that patients with high-PDE1B expression would have a better immune response to immunotherapies than those with low-PDE1B expression in osteosarcoma.
GSEA was applied by us to reveal PDE1B related pathways in osteosarcoma, and the outcomes of us showed that PDE1B was markedly linked to the calcium, cell cycle, chemokine, JAK STAT, and VEGF pathways.Previously articles showed that calcium homeostasis could inhibit cell proliferation and migration regulated by melatonin 32 .Zheng et al. identified POGZ as a hub gene in osteosarcoma, related to the cell cycle pathway 33 .Fan et al. found that dihydrotanshinone I could regulate osteosarcoma U-2 OS cells' adhesion and migration via chemokine pathways 34 .As for the JAK STAT pathway, Lv et al. utilized bioinformatics analyses and firstly  www.nature.com/scientificreports/mined serglycin as a prognostic marker involved with the JAK/STAT pathway 35 .Assi et al. summarized the roles of VEGF signaling pathways in osteosarcoma, pointing out future directions 36 .All of these five signaling pathways were significantly linked with osteosarcoma, and the PDE1B gene was predicted to play vital roles via these pathways in osteosarcoma.
To reveal the associations between PDE1B gene expression and immunity in osteosarcoma, immune checkpoints genes, tumor immune cells infiltration levels and tumor microenvironment were conducted based on the medium expression of the PDE1B gene, classified into high-and low-risk groups.As for immune checkpoints genes, our results reported that the PDE1B gene was significantly involved with immune checkpoints genes in osteosarcoma.Gul Mohammad et al. established an RBP-based model and found significant associations between this model and immune checkpoint genes, providing biomarkers to predict potential immune responses in osteosarcoma 37 .In the case of tumor immune cells infiltration levels, PDE1B was found to be markedly associated with T cells gamma delta and Macrophages M0 cells infiltration levels.As reported, Niu et al. found that Macrophages M0 and M2 cells were the major infiltrated cells in osteosarcoma by the CIBERSORT algorithm 38 .Wang et al. reported that the TUBB1 gene was positively linked to T cells gamma delta infiltration levels in osteosarcoma 39 .In terms of tumor microenvironment, our results presented that PDE1B gene expression was significantly associated with immune, ESTIMATE, and stromal scores by the ESTIMATE algorithm.Likewise in other published articles, Ma et al. successfully established an autophagy-related signature for risk stratification, survival prediction, and tumor microenvironment evaluation in osteosarcoma 40 .He et al. successfully constructed an immune-based LncRNA model to forecast osteosarcoma patients' prognosis and to evaluate the tumor microenvironment 41 .All of these indicated the vital associations between PDE1B gene expression and immunity in osteosarcoma.
However, no PDE1B targeted treatments were available currently.So, the TIDE algorithm was also applied by us to forecast immune responses by evaluating patients' PDE1B gene expression for immunotherapies for further clinical evaluations.The TIDE algorithm had been widely used in other studies for predicting immune responses.Wang et al. identified a ferroptosis-related model and applied the TIDE algorithm to show that lowrisk patients could benefit more from immunotherapies for lung squamous cell carcinoma 42 .Zhou et al. established a pyroptosis-based LncRNA model and also calculated the TIDE algorithm to forecast immune responses of their signature in kidney cancer 43 .In the present article, the outcomes of us shed light on that patients with high-PDE1B expression would have a lower TIDE score and a lower T cell dysfunction score.In other words, patients with high-PDE1B expression would have a better immune response to immunotherapies than those with low-PDE1B expression.All of these suggested that elevated PDE1B gene expression could prevent immune escape from osteosarcoma.
Future research would focus on using advanced computational models to identify novel mRNA or miRNA biomarkers for complex human diseases, which was crucial for identifying biomarkers in diseases like osteosarcoma.Recent studies had demonstrated that integrating various data sources and algorithms could significantly improve predictive accuracy [44][45][46] .For instance, the Rotation Forest for Essential MicroRNA Identification (RFEM), which combined multiple miRNA functional features, had shown notable improvements in miRNA identification accuracy 47 .These developments indicated that computational approaches utilizing data integration and multi-model strategies were vital for future mRNA or miRNA biomarker research.In our subsequent studies, we would further mine hub mRNA or miRNA biomarkers for osteosarcoma via various advanced computational models.
In this article, we firstly mined the prognostic values and immunological roles of the PDE1B gene in osteosarcoma with the assistance of experimental verifications, making our results much more reliable.This article provided some references for future PDE1B related research in osteosarcoma and was anticipated to provide operational targets for future clinically personalized treatment targets for osteosarcoma.The major limitation was that the current paper was only a preliminary exploration of the role of PDE1B in osteosarcoma, which was more based on theoretical analysis.Further in-depth in vitro and in vivo experiments in our subsequent articles were still required to verify these findings.

Conclusions
The PDE1B gene was found to be a tumor suppressor gene in osteosarcoma, related to OS prognosis and immunity.Univariate and multivariate Cox regression analyses indicated that the PDE1B gene had independent abilities in predicting OS for osteosarcoma.GSEA revealed that PDE1B was markedly linked to the calcium, cell cycle, chemokine, JAK STAT, and VEGF pathways.Further TIDE algorithms shed light on that patients with high-PDE1B expression would have a better immune response to immunotherapies than those with low-PDE1B expression, suggesting that elevated PDE1B gene expression could prevent immune escape from osteosarcoma.Further in-depth in vitro and in vivo experiments were still required to verify our findings.

Figure 1 .
Figure 1.The expression of PDE1B in osteosarcoma; (A) Venn diagram for gene filtration; (B) PDE1B gene expression in GSE28424 dataset; (C) PDE1B gene expression in GSE33382 dataset; (D) PDE1B gene survival analysis in TARGET osteosarcoma dataset; (E) qRT-PCR results of PDE1B gene in osteosarcoma cell lines; (F) Western blot results of PDE1B gene in osteosarcoma cell lines; (G) Bar chart of western blot results.

Figure 2 .
Figure 2. Boxplots of PDE1B gene expression distribution in (A) gender of female and male; (B) race of African, Asian, and White; (C) first event of relapse, censored or none; (D) disease at diagnosis of metastatic or nonmetastatic; (E) primary tumor site of arm/hand/pelvis, and leg/foot.

Figure 7 .
Figure 7. Correlations between PDE1B expression and immunity in osteosarcoma; (A) Proportions of tumor immune cells infiltration levels in tissues of TARGET osteosarcoma dataset; (B) 22 kinds of immune cells in high-PDE1B and low-PDE1B subgroups; (C) Correlation analyses between PDE1B and Macrophages M0 cells infiltration levels; (D) Correlation analyses between PDE1B and T cells gamma delta infiltration levels; (E) Tumor microenvironment evaluations; (F) Correlation analyses between PDE1B and ESTIMATE score; (G) Correlation analyses between PDE1B and immune score; (H) Correlation analyses between PDE1B and stromal score; *P < 0.05; **P < 0.01; ***P < 0.001; ns not significant.

Table 3 .
Gene sets enrichment analysis of PDE1B mRNA expression in Target osteosarcoma dataset.NES normalized enrichment score, NOM nominal, FDR false discovery rate.